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We discuss multi-dimensional generalizations of multicanonical algorithm, simulated tempering, 
and replica-exchange method. We generalize the original potential energy function Eq by adding 
any physical quantity V of interest as a new energy term with a coupling constant A. We then 
perform a multi-dimensional multicanonical simulation where a random walk in Eo and V space is 
realized. We can alternately perform a multi-dimensional simulated tempering simulation where a 
random walk in temperature T and parameter A is realized. The results of the multi-dimensional 
replica-exchange simulations can be used to determine the weight factors for these multi-dimensional 
multicanonical and simulated tempering simulations. 



Monte Carlo (MC) and molecular dynamics (MD) 
simulations of frustrated systems such as spinglass and 
biomolecular systems are very difficult because their 
free energy landscapes are rugged and long equilibra- 
tion time is necessary. In order to overcome this dif- 
ficulty, generalized- ensemble algorithms have often been 
employed (for reviews, see, for instance, Refs. [3, 0, [3, 0] ) • 
Generalized-ensemble algorithms are based on artificial, 
non-Boltzmann weight factors so that random walks in 
potential energy space and other variable space may be 
realized. Once an optimal weight factor is found, one 
makes a single long production run. From the results 
of this production run, one can reconstuct canonical, re- 
alistic ensembles for a wide range of temperature and 
other parameter values by the single-histogram [B| or 
multiple-histogram 0, @] reweighting techniques. Multi- 
wring 
are 

three of the most widely used generalized-ensemble algo- 
rithms. (ST is also referred to as the method of expanded 
ensemble llOU and REM is also referred to as parallel 
tempering [13].) I n this article, we present general for- 
mulations for multi-dimensional extensions of these three 
methods, where we generalize the original potential en- 
ergy function by adding any physical quantity of interest 
as a new energy term so that a random walk not only in 
the original potential energy space but also in the addi- 
tional energy space is realized. 

Let us consider the following generalized potential en- 
ergy function of a system in state x: 



muitipie-mstogram [y, yjj reweigntmg teenmques. iviv 
canonical algorithm (MUC A) 0, E| , simulated temper 
(ST) [53, El, and replica- exchange method (REM)pJ 



(1) 



Here, there arc L + 1 energy terms, E (x) and Vg(x) 
(I = 1, •••,£), and \W are the corresponding cou- 
pling constants for Vi(x) (we collectively write A = 
(A^ 1 ), • • • , A( L ))). The partition function of the system 



at fixed temperature T and A is then given by 

Z(T, A) = J dxexp{-/3E x (x)) 

= J dE dV 1 ■■■dV L n{E 0l Vx,---,V L ) exp (-0E X ) , 

(2) 

where [3 = l/k^T, k& is the Boltzmann constant, 
and n(Eo, Vi, ■ ■ ■ , Vl) is the multi-dimensional density of 
states. Here, the integral is replaced by a summation 
when x is discrete. 

The expression in Eq. ([1]) is often used in simulations. 
For instance, in simulations of spin systems, Eq(x) and 
Vi(x) (here, L = 1 and x = {Si, S2, ■ • •} stand for 
spins) can be respectively considered as the zero-field 
term and the magnetization term coupled with the exter- 
nal field A' 1 '. (For Ising model, E = -Jj2<ij> SiSj, 
Vi = — Si, and = h, i.e., external magnetic field.) 
In umbrella sampling [l4j | in molecular simulations, Eq(x) 
and Ve(x) can be taken as the original potential energy 
and the "biasing" umbrella potential energy, respectively, 
with the coupling parameter (here, x = {qi, • ■ ■ , q^} 
where qr i are the coordinate vectors of the i-th particle 
and N is the total number of particles). For the molecular 
simulations in the isobaric-isothcrmal ensemble, Eq(x) 
and Vi(x) (here, L — 1) can be respectively considered 
as the potential energy U and the volume V coupled with 
the pressure V. (Namely, we have x = {q 1 , ■ ■ ■ , q N , V}, 
E = U,Vi = V, and A^ = V, i.e., E x is the enthalpy 
without the kinetic energy contributions.) For simula- 
tions in the grand canonical ensemble with N particles, 
we have x = {q 1: ■■■ ,q N , N}, and E (x) and Vi(x) (here, 
L = 1) can be respectively considered as the potential en- 
ergy U and the total number of particles N coupled with 
the chemical potential [i. (Namely, we have Eq = U, 
Vi = N, and A^ = — /x.) We remark that generalized- 
ensemble algorithms in various ensembles are also dis- 
cussed in Refs. 15, 16j]. Moreover, we can introduce any 
physical quantity of interest (or its function) as the ad- 
ditional potential energy term Vg. For instance, Ve can 



2 



be an overlap with a reference configuration in spinglass 
systems, an end-to-end distance and a radius of gyration 
in molecular systems, etc. In such a case, we have to 
carefully choose the range of A*^ values so that the new 
energy term X^Vg will have roughly the same order of 
magnitude as the original energy term E . We want to 
perform a simulation where a random walk not only in 
the E space but also in the Vi space is realized. As 
shown below, this can be done by performing a multi- 
dimensional MUCA or ST simulation. 

We first describe the multi-dimensional MUCA sim- 
ulation which realizes a random walk in the (L + 1)- 
dimensional space of E (x) and Vg(x) (£ — I, • • • ,L). 
In the multi-dimensional MUCA ensemble, each 
state is weighted by the MUCA weight factor 
W mu (Eo, V\, ■ ■ ■ , Vl) so that a uniform energy distribu- 
tion of Eq, Vi, ■ ■ -, and Vl may be obtained: 

Pmu(E ,V 1 ,---,V L ) 

oc n(E , Vi,--, V L )W mu (E , V X ,---,V L ) = const , 

(3) 

where u(Eq, Vi, ■ • • , Vl) is the multi-dimensional density 
of states. From this equation, we obtain 



W mn (E a ,V lr --,V L ) oc 



1 



n{Eo,V x ,---,V L ) (4) 
= exp(~fj a E mu (E Q ,V u ---,V L )) , 

where in the second line we have introduced an arbi- 
trary reference temperature, T a — l/k^Pa, and wrote 
the weight factor in the Boltzmann-like form. Here, the 
" multicanonical potential energy 1 ' 1 is defined by 

E mu (E , Vi,---, V L ) = k B T a lnn{E , Vl, ■ • • , V L ) . (5) 

The multi-dimensional MUCA MC simulation can be 
performed with the following Metropolis transition prob- 
ability from state x with energy E^ = E + J2e=i X^Vg 
to state x' with energy E x ' = E ' + J2e=i A W U/ : 



w(x — > x') 



min 1 



1. 



W mu {E ',V 1 ',---,V L ') 

W mu (E Q ,V u ---,V L ) 
n{E Q ,V 1 ,---,V L ) • 

n(Eo',Vi',--,V L ') 



(G) 



An MD algorithm in the multi-dimensional MUCA en- 
semble also naturally follows from Eq. ([4]) , in which a reg- 
ular constant temperature MD simulation (with T = T a ) 
is performed by replacing the total potential energy E x 
by the multicanonical potential energy E mu in the New- 
ton's equations for the fc-th particle (k = 1, • • • , N) (see 
Refs. 17, 18l| for one-dimensional version): 



Pk 



dE mu (E ,V 1 ,---,V L ) 



(7) 



Secondly, we consider a multi-dimensional ST simula- 
tion which realizes a random walk both in temperature 



T and in parameters A. The parameter set A = (T, A) = 
(T, A^ 1 ), • • • , A^)) become dynamical variables and both 
the configuration and the parameter set are updated dur- 
ing the simulation with a weight factor: 



Wst(A) = exp (-/3£ A + /(A)) 



(8) 



where the function /(A) = f(T, A) is chosen so that the 
probability distribution of A is flat: 

P ST (A) <x 

J dEodV! ■■■dV L n(E , V lr --,V L ) exp (-f3E x + /(A)) 
= const . 

(9) 

This means that /(A) is the dimensionless ( "Helmholtz" ) 
free energy: 

exp (-/(A)) 

cx J dE dV! ---dVL n{E , V u ---,Vl) exp(-/3£ A ) . 

(10) 

In the numerical work we discretize the parameter 
set A in M(= M x Mi x • • • x M L ) different val- 



ues: A m = (T mo ,A m ) = (T mo , X ( m\, ■ ■ ■ , Xml), where 
mo = 1, • • • , Mo, mi = l,---,Mi (£ = 1, • • • , L). Without 
loss of generality we can order the parameters so that 



Ti < T 2 < ■ ■ ■ < T Mo and x\ e) < X^' < ■■■ < X 



An 



Mi 



(for 



each I = 1, • • • , L). The free energy / is now written as 

/m ,mi,"',mL — / C^mo ' A mi , • ■ ■ , X mL )■ 

Once the initial configuration and the initial parameter 
set are chosen, the multi-dimensional ST is realized by 
alternately performing the following two steps: 

1. A canonical MC or MD simulation at the 
fixed parameter set A m = (T mo ,A m ) = 



(T, 



A (1) • 



Ami ) is carried out for a certain 



steps with the weight factor exp(— [3 mQ E\ ). 

2. We update the parameter set A m to a new parame- 
ter set A m ±i in which one of the parameters in A m 
is changed to a neighboring value with the configu- 
ration and the other parameters fixed. The transi- 
tion probability of this parameter-updating process 
is given by the following Metropolis criterion: 

W (A m ->A m±1 ) =mm^ j (n) 

= min (1, exp (—A)) . 

Here, there are two possibilities for A m ±i, and we 
have A m± i = (T mo ±i, • • • , A™] , • • •) with 

A = (f3„ l(j ±i — Pmo) E^ — (/m ±l,mi,— ,itit ~ fmo,mi,---,m L ) , 

(12) 

for T-update, and A m ±i = (T mo ,- 
with 

,,, : | — X mf )Vl — {fm Q ,---,mi±l,--- ~ fm Q ,---,mi,---) , 

(13) 

for A^ -update (for one of I = 1, • • • , L). 



A = /3 mo (A 
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We remark that the random walk in Eq and in Vi for 
the MUCA simulation corresponds to that in j3 and in 
for the ST simulation: 



E <— >/3, 



(14) 



They are in conjugate relation. 

We can perform the multi-dimensional MUCA and ST 
simulations when we have optimal weight factors. How- 
ever, we do not know these MUCA and ST weight factors 
a priori and need to estimate them by short preliminary 
simulations. For one-dimensional version, three methods 
are well-known to obtain the weight factors: The first 
one is to use recursion formulas [2j, the second one is to 
use Wang-Landau methods and the third one is to 
use a short REM simulation and the multiple-histogram 
reweighting techniques [2(J 21, 22, 23|, 24 1. Here, we gen- 



eralize this third method to multi-dimensional versions 
(see also Refs. 

We use the multi- dimensional replica- exchange method 
(MREM) J2j| to determine the multi-dimensional MUCA 
and ST weight factors. The system for MREM consists of 
M non-interacting replicas of the original system in the 
"canonical ensemble" with M(= M x Mi x • • • x Ml) 
different parameter sets A m (m = 1,---,M). Because 
the replicas are non-interacting, the weight factor is given 
by the product of Boltzmann factor for each replica: 



W-MREM 



M 

n ex p (- 



(15) 



REM closely follows the ST procedures described 
above. In step 1, a "canonical" MC or MD simulation 
at the fixed parameter set is carried out for each replica 
simultaneously and independently for a certain MC or 
MD steps. In step 2, we exchange a pair of replicas i 
and j which are at the parameter sets A m and A m +i, 
respectively. The transition probability for this replica 
exchange process is given by 



w(A m <-> A m+ i) = min (1, exp(-A)) . 



(16) 



where we have 

A = {(3 mo - (3 mo+1 ) (E Xm (>1) - E Xm ( g M)) , (17) 
for T-exchange, and 

A = ^ mo (AW-A2 +1 )(F,(^)-^(^)) , (18) 

for A^ -exchange (for one of I = 1, •■■,L). Here, q'- 1 ' 
and stand for configuration variables for replicas i 
and j, respectively, before the replica exchange. Usually, 
Mo/2 or Mi/ 2 pairs of replicas corresponding to neigh- 
boring T or \W are simultaneously exchanged, and the 
pairing is alternated between the two possible choices, 



and (T 2 ,T 3 ),(T 4 ,T 5 ),--- 
and (A^AfMA^Af), 



i.e., (T 1 ,r 2 ),(r 3 ,r 4 ),- 

(AfU^MAW Af), 
respectively. 

To obtain the canonical distributions, the multiple- 
histogram reweighting techniques @, 0] are particularly 
useful. Suppose we have made a single run of the 
MREM simulation with M(= M x Mi x • • • x M L ) 
replicas that correspond to M different parameter sets 



A m (m = 1, • • ■ ,M). Let N„ 



,(E ,Vi,---,Vl) 



and n mo ,mi,--;m L be respectively the (L+ l)-dimensional 
potential-energy histogram and the total number of 
samples obtained for the m-th parameter set A m — 
(T mo , Xm\ , ■ ■ • , Xml ) . The multipje-histogram reweight- 
ing equations are then given by 

n(E ,Vi,---,V L ) 



E 



N„ 



XE ,V u ---,Vl) 



E 



L mo ,mi ,-■ ■ 



mo,mi,'"i m L 



exp ( f mo 

(19) 



and 



CX P( fm ,m 1 ,--,m L ) 

E n ( E o,Vi,---,VL)exp(-p mo E Xri 

E ,V t ,-,V L 

(20) 

The density of states n(Eo, Vi, ■ ■ ■ , Vl) and the dimen- 
sionless free energy fm ,m 1 ,---,m L are obtained by solv- 
ing Eqs. (fT9|) and ((20|) self-consistently by iteration. 
The canonical probability distribution at any temper- 
ature T = l/fcs/3 with any potential-energy param- 
eter value A is then given by P(E ,Vi, ■ ■ ■ ,Vl) = 
n(E ,V lr --,V L )exp(-PE x ). 

Finally, the weight factors for multi-dimensional 
MUCA (see Eq. (gj)) and multi-dimensional ST (see 
Eqs. © and pH]) ) are obtained from the generalized den- 
sity of states n(Eo, Vi, ■ ■ ■ , Vl) and the dimcnsionless free 
energy f mo , mi ,-, mL , respectively. 

As an example of the applications of the present formu- 
lations, we now present the results of a two-dimensional 
ST simulation. The system is a biomolecluar system 
studied in Ref. [2§| . We set E\ = Ep + AE'sol, where we 
have L = 1 in Eq. ([I]) and Eq = Ep is the conformational 

Sag P(Ep,Eso!) 
-6 
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FIG. 1: Canonical distributions P(Ep, Esol) with the 32 pos- 
sible parameter sets (T mo ,A mi ), which were obtained by the 
short two-dimensional REM simulation. 
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FIG. 2: Time series (a) and histogram J?(mo,mi) (b) of the 
parameter labels mo and mi for (T mo ,A mi ), which were ob- 
tained by the two-dimensional ST simulation. 

energy of the biomolecule and V\ = Esol is the solvent 
energy. The simulations were started from randomly gen- 
erated conformations. We prepared eight temperatures 
which are distributed exponentially between T\ = 300 
K and T 8 = 700 K (i.e., 300.00, 338.60, 382.17, 431.36, 
486.85, 549.49, 620.20, and 700.00 K) and four equally- 
spaced A values ranging from to 1 (i.e., Ai = 0, A2 = 
1/3, A3 = 2/3, and A4 = 1). The total number of replicas 
is then 32 (= 8 x 4). 

In Fig. 1, the canonical probability distributions at 32 
conditions obtained from the two-dimensional REM sim- 
ulation are shown. For an optimal performance of the 
REM simulation, there should be enough overlaps be- 
tween all pairs of neighboring distributions, which will 
lead to sufficiently uniform and large acceptance ratios of 
replica exchange. We see in Fig. 1 that there are indeed 
ample overlaps between the neighboring distributions. 

Using the results of this MREM simulation, we ob- 
tained the two-dimensional ST parameters / mo , mi (mo = 
1, • • • , 8; mi = 1, • • • , 4) by the multiple-histogram 
reweighting techniques (see Eqs. ijHJ), (fl"9|). and (|20|) ). and 
performed a two-dimensional ST simulation. 

The time series of labels of temperature T and param- 
eter A is shown in Fig. 2(a). The random walk in both 
T space and A space was indeed realized. The histogram 
of labels of T and A is shown in Fig. 2(b). We did get an 
expected flat histogram in T and A. 

Finally, we remark that once the weight factors for 
the multi-dimensional MUCA and ST are obtained, they 
can give the weight factors for lower-dimensional cases. 
For instance, the weight factor for the multimagnetical 
algorithm [261 ] can be obtained from that for the two- 
dimensional multicanonical-multimagnetical ensemble 
by integrating out the Eq variable (zero-field term). 
Likewise, the weight factors for multibaric-multithermal 
algorithm can be reduced to those for multibaric- 
isothermal ensemble and isobaric-multithermal ensemble 
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